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SUMMARY 


A  computer  program  written  in  FORTRAN  has  been  developed  to 
calculate  the  calibration  coefficients  of  a  six-component  strain  gauge 
balance.  The  method  of  Least  Squares  is  used  in  the  calculation  procedures 
which  do  not  require  the  components  of  the  balance  to  be  loaded 
independently.  The  procedures  are,  in  fact,  independent  of  the  way  the 
balance  is  loaded  as  long  as  the  supplied  load /output  data  provide  sufficient 
information  for  the  determination  of  the  coefficients.  The  program 
incorporates  several  statistical  estimations  fo^  the  evaluation  of  the 
performance  of  the  calibration  process.  It  can  therefore  be  used  as  an 
experimental  tool  for  the  investigation  of  different  loading  schemes  which 
may  be  used  to  calibrate  a  strain  gauge  balance.  ^ 
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a  matrix  derived  from  the  Least  Squares  criteria,  defined  by  equation  (7) 
coefficient  of  partial  correlation 
coefficient  of  calibration 
calibration  coefficient  matrices 

a  matrix  derived  from  the  Least  Squares  criteria,  defined  by  equation  (5) 

a  component  load  of  strain  gauge  balance 

vector  of  component  loads 

vector  of  “non-linear”  loading  terms 

total  number  of  sets  of  data 

degree  of  freedom  in  statistical  estimates 

index  of  summation 

a  component  of  strain  gauge  balance  output 
coefficient  of  multiple  correlation 
standard  deviation  of  Hi 
standard  error  of  estimates 
a  variable  in  regression  equation 
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1  Introduction 


In  wind  tunnel  operations,  calibrating  a  multi-component  strain  gauge  balance  is  a  process 
of  determining  the  coefficients  of  an  equation  which  describes  the  relationship  between  the 
applied  loads  and  the  strain  gauge  outputs  from  each  component.  This  involves  choosing 
an  appropriate  calibration  equation  which  best  describes  the  behaviour  of  the  particular 
balance.  The  balance  is  then  loaded  in  a  specific  manner  and  the  output  from  each  strain 
gauge  component  is  noted.  A  sufficient  number  of  loadings  lias  to  be  applied  to  the  system, 
depending  on  the  complexity  of  the  chosen  calibration  equation,  in  order  that  the  coefficients 
of  the  calibration  equation  may  be  totally  determined. 

Ideally,  a  perfect  strain  gauge  balance  would  have  each  strain  gauge  output  correspond¬ 
ing  to  one  and  only  one  load  component.  However,  in  a  real  situation,  it  is  impossible  to 
completely  eliminate  interactions  from  other  components.  Very  often,  significant  non-linear 
interactions  (functions  of  a  combination  of  load  components)  are  present  as  well  as  the  usual 
linear  ones  (functions  of  a  single  load  component).  For  a  well  designed  6-component  balance, 
it  is  usual  to  express  the  strain  gauge  output  of  the  ith  balance  component  (/?,-,  i  =  1, ...,  6)  as 
a  second  order  polynomial  function  of  all  the  components  of  the  applied  load  (II j,  j  =  1, ...,  6). 
Thus  the  form  of  the  calibration  equation  of  a  6-component  strain  gauge  balance  is 

II  i  =  +  Ci.iHi  +  ■  •  •  +  Ci,$Hs 

T  Ci.nHl  -f  Ci.22Ill  +  •  •  •  +  Ci.eelll 

+  CU2HlH2  +  Ci.13H1H3+-+Ci.s6HsH6,  i  =  1,...,6  (1) 

where  the  C’s  are  the  calibration  coefficients  determined  during  balance  calibration  and  Hi 
are  the  component  loads.  The  calibration  coefficients  of  each  component  i  may  be  classified 
as  follows: 

1.  “linear”,  e.g.  Ci.j  for  j  =  1 , ..., 6; 

2.  “load  squared”,  e.g.  Ci.jj  for  j  —  1,...,6; 

3.  “load  cross  product”,  e.g.  Ci.jk  for  j  =  1, .  .  . ,  5,  and  k  —  ( j  T  1 ),...,  6. 

However,  due  to  the  difference  in  characteristics  of  a  strain  gauge  in  tension  and  in  compres¬ 
sion,  an  odd  function  representation  is  more  appropriate  to  describe  its  behaviour  on  both 
sides  of  the  origin.  Hence  the  load/output  relationship  of  a  strain  gauge  balance  may  better 
be  accounted  for  by  the  inclusion  of  the  “load  cubed”  terms  in  the  calibration  equation.  Thus 
the  right  hand  side  of  equation  (1)  is  modified  with  the  addition  of  the  terms 

Ci.ll\II\  +  Ci. 222^2  f  . f  I'l.smll fi 

T  he*  set  of  six  calibration  equations  for  the  six  components  can  be  represented  in  matrix  form, 
with  the  calibration  matrix  partitioned  into  “linear”  and  “non-linear”  sub  matrices: 

(/<’!  --  \C\\\H |  i  K;2|[//-|  (2) 


I 


where  [/£]  is  the  vector  of  strain  gauge  balance  component  outputs, 

[Cl]  is  the  matrix  of  linear  calibration  coefficients,  Ci.j, 

\C 2]  is  the  matrix  of  non-linear  calibration  coefficients,  C;  jj,  Cijk  and  Ci.jjj, 

[//]  is  the  vector  of  component  loads,  and 

[//']  is  the  vector  of  squares,  cross  products  .and  cubes  of  component  loads. 

The  conventional  method  of  calibrating  a  strain  gauge  balance  of  the  type  given  by 
equation  (2)  is  to  apply  pure  component  loads  independently.  When  a  single  component 
is  loaded,  the  outputs  of  all  components  of  the  strain  gauge  balance  respond  to  one  and 
only  one  load  component.  The  “linear”,  “load  squared”  and  the  “load  cubed”  coefficients 
may  hence  be  determined.  The  “load  cross  product”  coefficients  are  found  by  applying  two 
pure  component  loads  simultaneously  to  the  balance  keeping  one  load  constant  at  a  finite 
value  while  varying  the  other.  This  method,  although  straightforward  to  use,  relies  on  pure 
component  loads  which  rarely  represent  the  real,  practical  loading  situation  which  the  model 
encounters.  It  contravenes  the  basic  principle  that  the  loadings  of  the  balance  during  the 
calibration  should  be  representative  of  that  experienced  by  the  model  during  wind  tunnel 
tests. 

Theoretically,  if  the  calibration  equation  describes  exactly  the  relationship  between  the 
applied  load  and  the  strain  gauge  output,  the  way  in  which  the  balance  is  loaded  during 
calibration  should  not  affect  the  outcome  of  the  results.  Whether  the  balance  components 
are  loaded  independently  or  simultaneously,  the  values  of  the  calibration  coefficients  should 
be  the  same.  However,  in  reality,  the  calibration  equation  given  by  equation  (2)  is  only  an 
approximation  and  does  not  describe  exactly  the  load/output  relationship  of  a  practical  strain 
gauge  balance.  The  responses  of  the  strain  gauge  balance  components  can  be  quite  different 
under  different  loading  distributions  which  causes  different  deflection  characteristics  in  the 
balance.  The  calibration  coefficients  thus  calculated  are  dependent  of  how  the  balance  is 
loaded  during  the  calibration. 

Ramaswamy  et  al.  [1]  proposed  a  Least  Squares  (LS)  method  for  the  calibration  of 
wind  tunnel  balances.  This  method  does  not  require  the  balance  components  to  be  loaded 
independently.  In  fact,  any  combination  of  component  loads  may  be  applied  simultaneously. 
The  distributions  of  loads  on  the  balance  during  the  calibration  can  be  chosen  to  match  as 
closely  as  possible  those  likely  to  be  experienced  by  the  model  during  wind  tunnel  tests.  In 
addition,  there  is  no  restriction  of  the  sequence  or  pattern  in  which  the  loads  have  to  be 
applied  in  order  to  obtain  a  particular  calibration  coefficient  .  Virtually  any  random,  set  of 
component  loads  may  be  applied  to  the  balance  in  any  sequence  provided  that  all  the  data 
so  obtained  generate  enough  independent  information  for  the  calculations  of  the  coefficients. 
Another  advantage  of  the  LS  method  is  that  the  “linear”  and  “non-linear”  coefficients  can  be 
determined  at  the  same  time.  This  is  because  by  loading  all  the  components  simultaneously, 
the  “linear”  and  “non-linear”  interactions  among  the  components  are  all  accounted  for  in  the 
data.  Less  effort  will  hence  be  required  to  calibrate  the  balance  if  one  can  design  an  optimum 
loading  scheme  which  would  simulate  the  real  loading  situation  as  well  as  providing  the  most 


information  for  the  determination  of  the  coefficients  with  the  least  number  of  sets  of  loads. 


A  major  shortcoming  of  the  LS  method  is  that  it  assumes  the  independent  variables,  i.c. 
the  IPs  in  equation  (1),  are  free  of  noises.  In  reality,  due  to  certain  factors  such  as  friction 
or  misalignment  of  the  point  of  application  of  //,,  even  though  the  same  load  is  applied  to 
the  balance  twice,  the  actual  load  experienced  bv  the  balance  may  not  be  exactly  the  same 
in  both  cases.  Nevertheless,  the  errors  arising  from  the  difference  between  the  load  applied 
and  the  actual  load  experienced  by  the  balance  are  usually  insignificant  when  compared  with 
the  noises  in  the  outputs  R{.  Furthermore,  for  a  well  designed  strain  gauge  balance,  the 
noises  in  R{  are  often  satisfactorily  small  such  that  the  LS  method  can  normally  be  applied 
successfully. 

This  report  describes  a  computer  program,  SGBCalib,  which  implements  the  LS  method 
of  determining  the  coefficients  of  the  calibration  equation  given  by  (2).  The  only  inputs  to  the 
program  are  the  sets  of  applied  component  loads  and  the  corresponding  strain  gauge  outputs. 
If  the  given  data  contain  sufficient  information,  the  calibration  coefficients  will  be  evaluated. 
The  program  also  generates  useful  statistical  estimates  of  how  closely  the  resulting  calibration 
equation  fits  the  given  data.  The  reliability  of  the  calibration  equation  is  expressed  in  terms 
of  the  standard  error  of  estimate  of  the  strain  gauge  outputs  for  all  the  six  components. 
The  95%  confidence  levels  of  the  “linear”  coefficients  are  also  estimated  as  indications  of  the 
accuracy  of  the  calibration. 

Apart  from  the  function  of  calculating  the  calibration  coefficients  from  a  given  set  of 
data,  SGBCalib  can  be  used  to  assist  in  the  investigation  of  how  the  balance  components 
should  be  loaded  so  as  to  yield  the  most  information  with  the  least  number  of  sets  of  applied 
loads  and  at  the  same  time  best  simulate  the  patterns  of  loading  expected  in  the  wind  tunnel 
testings.  It  can  also  be  used  for  studies  of  whether  the  errors  of  a  particular  strain  gauge 
output  could  be  reduced  by  further  loadings  designed  to  optimally  minimize  the  uncertainties 
in  the  calibration  coefficient. 


1  Formulations  of  the  Least  Squares  Method 


Given  that  for  a  set  of  component  loads,  Hi  where  i  =  1,  ...,6,  the  corresponding  strain  gauge 
balance  outputs  are  R.{.  IF  there  are  A'  sets  of  //,  with  N  sets  of  corresponding  Rt)  denoted 
by  IIt,p  and  /£;  p  where  p  —  1, ...,  N ,  then  the  calibration  coefficients  may  be  found  such  that 
the  sum  of  the  squares  of  the  difference  between  the  measured  strain  gauge  output  and  that 
obtained  from  the  calibration  equation  is  a  minimum.  That  is 

N  , 

f  i  —  1 1'  i  1  fil.p  i  fi.2  IL.p  1  1  f  i  GGfi  H  G'jj  ' '  R  i  ^  ■  1  ,  ■  -  ■ ,  5  , 

P  ;  • 

should  be  a  minimum.  This  involves  partial  different iating  each  of  the  >:,  with  respect  to  each 
of  tin’  coefficient  s  and  equaling  the  resulting  expressions  to  zero.  The  result  is  a  set  of  TS 


non-linear  simultaneous  equations  for  each  component  i  given  as  follows1: 


y  ,  [Ci.,//,  4  Ci.jII 2  4  ■ 

•  f  C, .66g//g  “  «;]  //> 

-  0 

Y,  [Co/fi  +  C{  2  U  2  4  • 

■  ■  +  Ci60G//3  -  Ri]  Il2 

=  0 

[C,.,//,  4  6Y2//2  + 

4-  C,.GGG//3  -  Rt ]  III 

=  0. 

(3) 

It  is  more  convenient  to  put  equations  (3)  in  matrix  form: 

[£)[C,]  =  \Ai)  ,  *=  (4) 


where 


(£]  = 


E«2//i  E  W2//2 


E^.^6 

E  ^2//c3 


E^l  E'^2  •  E  Ul  "l  J 


(5) 


is  a  33  X  33  square  matrix  and  is  common  for  all  the  six  components.  The  two  column 
matrices,  [C;]  and  [/!,],  for  the  ith  component  are  given  by 


~  Cm 

'  E^i  Ri  ' 

[c,]  = 

Ci.2 

Mi]  = 

E  H2  Ri 

Ci.  66C) 

.  E  H\ }r{ 

The  [Cfj’s  and  [d,)’s  for  all  the  six  components  may  be  combined  to  form  two  33  X  6  matrices, 
[C]  and  [/!],  such  that 


[Cl 


Ml 


Cm 

C2.1 

CG.i 

C,.2 

C2.2 

CG.  2 

C 1  GGG 

C2  GGG 

f'CTiGr, 

E//i« 

1  M 

//. 

//, 

Hr, 

E  h  2  R 

.  E 

//-, 

R  2 

V' 

z_> 

Hz 

/fc, 

E  //./«• 

,  E 

R-, 

V 

Hi 

/Ai 

(6) 


(') 


'For  the  sake  of  clarity,  the  subscript  p  for  It,  ami  II,  is  omitted  in  the  remainder  of  the  report.  Whenever 
the  summation  symbol  is  used,  it  is  understood  that  the  sum  is  to  be  taken  over  p  1,  —  A 


(8) 


Tlie  Least  Squares  criteria  can  therefore  he  represented  by  tiie  matrix  equation 

m\c\  mi- 

Provided  that  | K )  is  non-singular,  the  coefiicient  matrix  [(7)  may  he  found  as 

K-’i  -  I/--T ‘Ml  (fl) 

3  Statistical  Estimations 


In  this  section  we  adopt  a  few  statistics  to  measure  the  accuracy  and  reliability  of  the  cal¬ 
ibration  of  the  balance  using  the  LS  method.  We  seek  to  determine  how  well  the  resulting 
calibrating  equation  describes  the  relationship  between  the  strain  gauge  outputs  and  the  ap¬ 
plied  component  loads.  The  statistics  described  below  provide  us  with  information  on  and 
measurements  of  the  goodness  of  fit  of  the  calibration  equation,  and  the  uncertainties  of  the 
coefficients  determined  by  the  LS  method. 


3.1  Standard  Error  of  Estimates 


Tli<?  standard  error  of  estimates  of  the  ith  component  is  given  by 


Sitn,i 


E(fl.  -  RiY 

N, 


(10) 


wher.e  /t,  is  the  measured  ith  component  of  the  strain  gauge  balance  output. 

R,  is  the  corresponding  estimated  value  using  the  calibration  equation. 

Nf  is  the  degree  of  freedom,  which  is  (fV  -  33),  and 
N  is  the  total  number  of  sets  of  loading  cases. 

It  is  a  measure  of  the  scatter  of  the  balance  outputs  about  the  estimated  values.  It  has 
properties  analogous  to  those  of  the  standard  deviation,  i.e.,  if  N  is  large  enough,  there 
would  be  approximately  68%,  95%  and  99. 7%  of  the  measured  Rt  lies  within  the  bounds  of 
2.S’/f//,  and  3.S' respectively  from  the  estimated  R,  value. 


3.2  Coefficient,  of  Multiple  Correlation 

The  coefficient  of  multiple  correlation  of  the  ft  h  component,  is  defined  as 


' 


>;(/c  A’,)2 

v(/>\  /<\)2 1  R,)' 


The  term  V( R,  II,)2  is  known  as  the  explained  variation ,  so  called  because  it  accounts  for 
the  variations  (of  R,  from  the  mean)  which  are  explained  by  the  regression.  £](/?,  -  R{)2 
is  the  unexplained  variation.  It  represents  the  variations  which  remain  unexplained  after 
the  regression.  The  coefficient  of  multiple  correlation  shows  what  fraction  of  the  overall 
variations  is  accounted  for  by  applying  the  calibration  equation.  It  is  therefore  a  measure 
of  the  correlation  between  a  balance  component  output  and  all  of  the  component  loads.  It 
shows  how  well  the  calibration  equation  describes  the  relationship  between  the  outputs  of 
the  strain  gauge  balance  and  the  component  loads.  It  is  also  a  measure  of  the  closeness  of  fit 
of  the  calibration  equation  to  the  actual  measured  strain  gauge  balance  outputs.  The  value 
of  r,  lies  between  0  and  1.  The  closer  it  is  to  1  the  better  is  the  linear  relationship.  When 
r,  1  the  correlation  is  called  perfect. 


3.3  Coefficient  of  Partial  Correlation 


'Phe  coefficient  of  partial  correlation  measures  the  linear  relationship  between  a  component 
of  the  strain  gauge  balance  output  and  a  single  component  load  while  the  other  component 
loads  are  kept  constant. 


The  linear  correlation  coefficients  a\2  between  two  variables  A’|  and  A’ 2  is  given  by 

flj2  ^-^-(EA'^A-;) _  (12) 

vfa  £  A',2  -  ( E  A' J  )2 1  ( jv  £  A'f  -  ( E  A  2 )"2 1 

where  A'  is  the  number  of  data  points.  a!2  is  called  the  zero  order  correlation  coefficient.  We 
can  denote  the  coefficient  of  partial  correlation  between  A’i  and  A'2  keeping  A’ 3  constant  by 
a 1 2  3  (the  1st  order  correlation  coefficient)  such  that 


a12  a  1 3fl23 

'll  2.3  -  -  r-  ---  - 

V  (  '  “  *1 1  3  )  (  *  a2.s) 


(13) 


Similarly  the  coefficient  of  partial  correlation  between  A’i  and  A2  keeping  A’ 3  and  A’ 4  constant 
may  be  given  by  (see  fir  e  g.  Spiegel  [2]) 


"12-1  "l3.-1rt23  4  "12.3  "M3"21.3 

"12.31  j  - - '  ,  -  -  -  (Id) 

\/(i  «^.,)(i  "23.1)  \/  ( 1  "ii.iH1  "5.1.3) 

The  coefficient  of  partial  correlation  between  the  strain  gauge  balance  output  R,  and  the 
component  load  // 1 ,  sav,  keeping  the  other  component  loads  constant  may  therefore  be  de¬ 
noted  by  (?/;  i  224 r,r,,  where  the  numerical  subscripts  denote  the  component  loads,  and  can  be 
deduced  from  equal  ions  ( 12),  ( I  3)  and  (It) 


An  ideal  strain  gauge  balance  would  have  perfect  linear  correlation  (u  I)  between 
an  output  and  its  corresponding  component  load  and  zero  correlation  with  the  rest  of  the 
components  These  coefficients  therefore  indicate  how  lar  the  strain  gauge  balance  has  fallen 
short  of  t  Ins  ideal  behaviour 


(1 


3.4  95%  Confidence  Interval 


To  establish  a  95%  confidence  interval  for  a  calibration  coefficient  C \.j  we  can  construct  the 
t -distribution  statistics: 


<  = 


Ci.j  -  Cj.j 

Snn,t/SH'j 


i  "  j  =  I,..-, 6 


where  Sjih, i  is  the  standard  error  of  estimate  of  the  ith  component  given  in  section  3.1, 
Sn,j  is  the  standard  deviation  of  II  j, 

Ci.j  is  the  expected  mean  of  Ci.j. 

The  95%  confidence  interval  for  the  coefficient  G\  j  is  then  given  by  the  expression 


<0.975  (  S [ill ,i  \ 

Nf  \  Sn.i  J 

The  percentile  values  <0.975  f°r  various  degree  of  freedom  (N  f)  may  be  found  in  most  reference 
books  in  statistics  (see  e.g.  Yamane  [3]). 


4  The  Computer  Program  SGBCalib 

The  program  SGBCalib  was  written  in  FORTRAN  and  was  implemented  in  DEC1  MicroVAX1 
II  operating  under  the  VMS1  (Version  4.5)  system.  The  program,  however,  was  not  intended 
to  be  VMS-specific,  and  the  language  used  conforms  to  the  American  National  Standard 
FORTRAN-77  (ANSI  X3.9-1978). 


4.1  Program  Variables 

Drag  measurements  using  wind  tunnel  force  balance  usually  require  high  order  of  accuracy.  In 
order  to  maintain  data  integrity,  all  the  program  variables  concerned  are  assigned  as  “double 
precision”.  Most  of  the  variables  were  passed  between  the  various  subroutines  via  common 
blocks: 

/matrixA/,  /matrixC/,  /matrixE/,  /inv/ 

These  common  blocks  contain  the  elements  of  the  matrices  [.4],  [C], 

| /•;)  and  | A’ |  '*  respectively. 

/loads/, /readings/ 

1  I)K(\  MicroVAX  nnH  VMS  ;»rc  registered  trademarks  of  Digital  Dtpument  < 'nrpnmtion 


They  store  all  the  given  //,  and  //,  data  in  the  2-D  arrays  HH  and  RR 
respectively.  The  program  at  present  allows  a  maximum  of  200  sets  of 
//,  and  R.{  values  to  be  processed.  This  may  be  changed  if  necessary 
by  altering  the  dimension  of  the  second  subscripts  of  the  statement 

double  precision  HH(6,200),  RR(6,200). 

The  dimension  of  the  second  subscript  of  the  variable  res (6, 200) 
in  the  subroutines  Statistics,  SumStats  and  Chart  should  also  be 
changed  to  be  the  same  as  those  of  HH  and  RR. 

/Hvector/  This  is  an  array  of  the  33  //  terms  in  the  calibration  equation.  That 
is,  it  contains  the  values  of 


//1,//2,...,//6,//12,//22, 

//  3  //  3  H  3 

11  1  i  n  2  i  '  •  •  i  "6  • 


.,//g2,//1//2,//1//3,...,//5//6, 


/ residual/ 


/misc/ 

/ statsvar/ 


This  common  block  appears  in  the  statistical  routines  only.  It  consists 
of  the  residual  values  of  Rit  i.e.  ( Ri  -  Ri).  The  block  also  contains 
the  maximum  value  of  the  residual,  resmax,  which  is  used  for  scaling 
of  the  chart  of  the  residuals. 

This  is  a  miscellaneous  block  containing  the  pth  set  of  values  of  Ri,  R, 
and  Hi.  The  mean  values  of  Ri  are  also  included  in  this  block. 

This  block  contains  the  variables  used  for  the  calculations  of  the  vari¬ 
ous  statistical  estimates  for  the  calibration.  The  major  variables  are: 


sRH(i , j ) 

=  HRiHj, 

sR(i) 

=  ZRi 

sR2(i) 

= 

sH(i) 

=  Y.IU 

sH2(i) 

=  £  //* 

sHH(i, j) 

=  ZHiHi, 

dRRe2(i) 

=  £(«<-«.-)' 

dRRm2(i) 

dReRra2(i) 

=  E{Ri-Ri) 

4.2  Program  Structure 

The  program  is  made  up  of  different  subroutines  designated  for  various  specific  tasks.  The 
listing  of  these  subroutines  is  given  in  Section  .r>.  The  tasks  undertaken  by  these  subroutines 
may  be  divided  broadly  into  three  areas: 

1.  Program  inputs  and  outputs. 


8 


2.  Calculations  of  the  matrices  [A],  [E j,  [/?]  1  and  \C\. 

3.  Estimation  of  the  statistics  of  the  calibration. 


Apart  from  the  various  subroutines  written  specifically  for  the  LS  method  of  calculating  the 
calibration  coefficients,  there  are  two  main  subroutines  which  are  used  in  matrix  operations. 
They  are  the  matrix  inversion  (Matlnv)  and  the  matrix  multiplication  (MatMult)  routines. 
For  the  sake  of  completion,  these  subroutines  are  also  listed  in  Section  5. 


4.2.1  Program  Inputs  and  Outputs 

The  only  program  inputs  are  the  sets  of  Hl  and  R{  values  which  are  read  from  the  data  file 
SGBIN.DAT.  The  format  of  the  file  is  as  follows: 

1.  The  first  line  is  empty  and  is  skipped  by  the  program.  It  may  be  used  for  data  identi¬ 
fication  purpose. 

2.  The  second  line  contains  a  row  of  6  components  of  loads,  H{,  each  of  which  is  separated 
from  the  next  value  by  one  or  more  spaces. 

3.  The  third  line  is  a  row  of  6  components  of  strain  gauge  outputs,  Ri,  corresponding  to 
each  of  the  //;  value  in  the  above  row. 

4.  This  pattern:  an  empty  line,  a  row  of  6  /f,  values  followed  by  a  row  of  6  corresponding 
Ri  values,  repeats  throughout  the  data  file  for  as  many  sets  of  loading  as  are  available. 

When  the  program  is  executed,  all  the  relevant  information  and  results  are  written  o 
an  output  file  SGBOUT.DAT  which  may  subsequently  be  printed  on  a  line  printer.  The  outpu  : 
of  the  program  consist  of  the  following: 

1.  All  the  values  of  //,  and  Rl  provided  in  the  input  file. 

2.  The  evaluated  calibration  coefficient  matrix  [C|  in  the  form: 


Ci.i 

C\  2 

C  2.1  C-S  1 

C2.2 

C4.1 

C5.1 

c  0.1 

E*o  2 

Ox.u 

Ci.u 

Eo.i  1 

(■  1  50 

<■■2.50 

Efl.so 

('  1  000 

E  2.000 

( -0.000 

3.  Statistical  estimations  which  include: 

•  the  standard  errors  of  estimates  for  each  of  the  component, 

•  the  coefficients  of  multiple  correlation, 

•  the  coefficients  of  partial  correlation  between  and  each  of  the  component  loads, 

•  the  95%  confidence  level  of  the  linear  coefficients,  i.e. 


Cm 

Cj. 2 

C2.1 

C2.2 

C3.1 

C4.1 

Cs. 1 

Cr,i 

Cg.2 

Cl. 6 

C2.6 

C3.G 

C4.6 

C5.6 

C6.6 

4.  A  table  of  comparisons  of  the  measured  and  estimated  values  of  Ri  and  the  residuals. 

5.  A  line  printer  plot  of  the  residuals  of  Ri  over  Ri  for  each  of  the  6  components.  The 
maximum  deviation  is  ±10  characters  from  the  zero  marker  (the  character  “!”)  and  is 
scaled  to  be  equivalent  to  the  maximum  residual  value  of  all  6  components. 


4.2.2  Evaluations  of  the  Matrices 

As  each  set  of  the  Ri  and  Ri  values  are  read  from  the  input  file,  the  “load  squares”  (//?), 
“load  cubes”  (Hf)  and  “load  cross  products”  ( HiHj ,  i  =  1,.  ..,5  and  j  -  i  +  1,. .  .,6)  are 
evaluated  by  means  of  the  subroutine  Hf  actors.  The  elements  forming  the  matrices  \E\  and 
[A]  are  then  calculated  and  summed  in  the  subroutine  SumElements. 

After  all  the  Hi  and  Ri  values  have  been  read,  [£]_1  is  found  using  the  matrix  inversion 
routine  Matlnv.  An  error  flag  (designated  by  the  variable  error)  is  used  to  indicate  if  an 
inverse  has  successfully  been  calculated.  When  error  —  0,  the  inverse  is  found  and  no  error 
occurs.  When  error  =  2,  \E)  does  not  have  an  inverse  and  the  program  terminates  with  the 
message 

INVERSE  OF  [E]  DOES  NOT  EXIST  HI 

appearing  on  the  terminal  screen.  Wlten  [/J]-1  is  found  tlu>  coefficient  matrix  [C|  is  evaluated 
according  to  equation  (9)  by  the  matrix  multiplication  routine  MatMult  and  the  elements  of 
\C\  are  written  to  the  output  file  SGB0UT.DAT. 


4.2.3  Statistical  Estimations 

The  evaluations  of  the  various  statistics  ment  ioned  in  Sect  ion  3  are  initialed  by  t  lie  sub  rout  ine 
Statistics.  The  estimated  value  of  denoted  as  Ro(i)  in  the  program,  is  calculated 


It) 


according  to  equation  (2).  Before  equation  (2)  can  he  applied,  the  coefficient  matrix  [C]  is 
transposed  and  partitioned  into  a  linear  matrix  [Cl]  and  a  non-linear  matrix  [C2j.  This  is 
accomplished  by  the  statements 

do  i  =  1 ,  6 
do  j  =  1 ,  6 

c 1 (i ,  j )  =  c(j ,i) 
end  do 

do  j  =  1,  27 

c2(i  ,  j )  =  c(j+6,i) 
end  do 
end  do 

in  the  Statistics  routine.  The  subroutine  SGBReading  then  evaluates  Re(i)  from  the  data 
of  Hi  according  to  equation  (2)  and  the  various  summations  required  by  the  statistics  are 
calculated  in  the  subroutine  SumStats. 

The  various  statistics  described  in  Section  3  are  designated  as 

S(i)  for  the  standard  error  of  estimates  5/?//,;, 

cmc(i)  for  the  coefficients  of  multiple  correlation  t\,  and 

cl(i,j)  for  the  95%  confidence  levels  of  the  linear  coefficients  Cij. 

Because  of  the  complexity  of  the  calculations  of  the  coefficients  of  partial  correlation,  they 
are  described  in  more  detail  in  Section  4.2.4.  The  calculation  of  the  confidence  levels  requires 
the  0.975  percentile  values  at  various  degree  of  freedom  for  the  t-statistics.  These  values  are 
provided  as  a  look-up  table  in  the  subroutine  Tstats. 

A  graphical  representation  of  the  residuals  of  Hi  over  /<!,•  is  constructed  by  the  subroutine 
Chart.  The  results  for  all  six  components  are  displayed  together  across  the  width  of  the  page. 
The  zero  axes  are  represented  by  the  character  “!”.  The  residuals  for  components  1  and  4 
are  represented  by  the  character  those  for  components  2  and  5  by  the  character  “o"  and 
those  for  components  3  and  6  by  the  character  “#”•  The  scale  of  the  deviation  is  based  on 
the  maximum  value  of  the  residual  (designated  by  the  variable  resmax  in  the  program)  of  all 
six  components.  The  value  of  resmax  occupies  a  space  of  10  character  on  the  chart  and  the 
other  values  of  the  residual  are  scaled  accordingly. 

4.2.4  Evaluation  of  Coefficient  of  Partial  Correlations 

For  t  he  sake  of  clarity,  let  us  consider  a  single  component  of  the  strain  gauge  balance  output 
ll:\,  say,  and  seek  to  determine  the  correlat  ion  coefficients  between  and  each  of  t  he  compo¬ 
nent  loads  //,  keeping  the  other  component  loads  constant.  It  is  convenient,  then,  to  consider 


an  array  A',,  i  =  such  that  A'i  -  /i3,  X2  =  // 1,  A3  =  //2 ,  ...,  and  A'7  =  //G.  The 

coefficients  of  partial  correlation  between  Ii2  and  each  of  the  //’ s  may  therefore  be  denoted 
as 

«12. 34567,  «13. 24567,  <<14 .23567,  <115.23467,  <<16.23457,  aild  <1 17.23456- 

To  determine,  say,  <112.34507,  we  may  use  the  relationship 


«12. 34567  = 


<112.4567  ~  <ll3.4567<l23.4567 


A 


1 


13  4567 


)(1- 


<23.4567 


(15) 


in  which  the  variable  A'3  is  used  as  the  “mediator”.  Theoretically,  we  can  use  any  of  the 
“kept-constant”  variables  A'3,  A4,  A'5,  A'e  and  A'7  to  be  the  mediator  to  calculate  <112. 34567- 
For  convenience  in  programming,  we  opt  for  the  lowest  subscript  of  the  “kept-constant” 
variables. 

To  determine  the  next  lower  order  correlation  coefficient,  say,  <113.4507  in  the  right  hand 
side  of  equation  (15),  we,  therefore  use  the  relationship 


<113.4567 


<113.567  ~  <tl4.567<l34.567 
\/(*  ~  a14.S67)(l  ~  a34.567) 


(16) 


in  which  is  now  the  mediator.  Ultimately,  these  coefficients  are  related  to  the  zero  order 
correlation  coefficients  an,  ai3,  ...,  a2 3,  a24,  ...,  a67  where 


NY  AAA',  -  (E  A'.XE  Xj) 
y/l*± A? -(^i)2PE  Xj  -  (E  Xj)2 


(17) 


In  the  program,  the  determination  of  the  coefficients  of  partial  correlation  is  carried  out 
by  the  subroutine  PartCorr.  The  sums  Y  A',,  Y  A'f,  and  A UXj  are  mapped  into  the  arrays 
suinX(i)  ,  sumX2(i),  and  sumXX(i)  respectively.  The  correlation  coefficients  are  denoted  by 
the  3-D  array  corr(m,i,  j)  in  which  (m-1)  is  the  order  of  the  correlation,  i  is  the  subscript 
of  the  dependent  variable  A{  and  j  is  that  of  the  independent  variable  Xj.  The  values  of  the 
correlation  that  we  are  looking  for  are  therefore  stored  in  corr(6,l,k)  ,  k  =  2,  ...,7. 

The  zero  order  correlation  coefficients  are  evaluated  first  and  are  determined  by  the 
statements 


do  ml  =  1,6 

do  m2  =  (ml  +  1 )  , 7 

corr ( 1 ,ml ,m2)  =  (N*sumXX (ml ,m2) -sumX (ml ) *sum.X (m2) ) 

*  /sqrt  (  (N*sumX2(ml. )  -  sumX  (ml)*sumX(ml))* 

*  (N*5umX2(m2)-5umX(m2) *sumX(m2) ) ) 


end  do 


The  mediator  for  the  determination  of  each  order  of  the  correlation  coefficients  is  pre¬ 
selected  in  the  subroutine  Rank.  The  “ranking”  of  the  mediator  subscript  is  stored  in  the 
array  mediator  (m)  so  that  the  highest  order  (i.e.  m  =  6)  would  have  the  lowest  subscript  in 
the  range  of  2  to  7  excluding  the  value  of  k  when  determining  corr  (6 , 1 ,  k) .  The  next  highest 
order  would  then  have  the  next  lowest  subscript  and  so  on.  For  example,  to  determine  the 
coefficient  aj5 .23467)  the  "ranking”  of  the  mediator  subscripts  is 

order  :  2  3  4  5  6 

mediator  :  76432 

To  calculate  the  1st  order  correlation  coefficients,  the  variable  with  the  subscript  7  is  to  be 
the  mediator.  Hence  the  coefficients  an. 7,  <113.7,  <114.7,  •  .  . ,  <156.7  are  evaluated. 

Because  of  the  symmetrical  relationship  that  corr(m,i,j)  =  corr(m,j,i),  only  the 
values  of  corr(m,i,  j)  in  which  i  <  j  are  evaluated  and  used  in  the  calculations. 

5  Program  Listing 

The  following  pages  are  a  full  listing  of  the  program  SGBCalib.  The  program  is  intended  to 
be  executed  interactively  on  a  computer  terminal.  It  may  also  be  submitted  as  batch  job  in 
most  system  without  any  modification. 

The  program  listing  is  roughly  divided  into  3  sections: 

Main  section  This  section  consists  of  the  main  routine  SGBCalib  and  the  subroutines 
SumElements  and  Hf actors.  They  are  the  routines  that  implement  the  LS  method. 

Statistical  section  This  includes  the  subroutines  Statistics,  Statslnitialize,  SumStats, 
SGBReading,  PartCorr,  Rank,  Chart  and  Resetline  and  are  responsible  for  the  evalu¬ 
ations  of  the  statistical  variables. 

Matrix  Operations  The  matrix  inversion  Matlnv  and  matrix  multiplication  MatMult  rou¬ 
tines  are  included  for  the  sake  of  completion  of  the  program  listing.  There  are  other 
minor  subroutines  following  Matlnv  and  are  called  by  Matlnv  only. 


Program  SGBCalib 


Strain  Gauge  Balance  Calibration  using  Least  Squares  Method 

Given  a  set  of  S.G.B.  readings,  R(i),  from  a  set  of 
known  loadings,  H(i),  the  matrices  [E]  and  [A]  are  formed. 
The  coefficient  matrix,  [C] ,  is  obtained  according  to 
the  equation 

[C]  =  [Einv]  [A] 

Written  by  S.  Lam  -  31  March  1988 

Associated  routines:-  Hfactors(H) 

SumE 1 ement  s ( R ) 

Statistics (N) 

Hf actors (H) 

Matlnv(dim,  E,  Einv,  error) 

Mat Mult (Cl ,m,n,C2,n,k,C3, error) 


double  precision  E(33 ,33) ,Einv(33 ,33) , A(33 ,6) ,C(33 ,6) 
double  precision  H(6) ,R(6) ,Hstar(27) ,Rm(6) 
double  precision  sumR(6) ,HH(6 ,200) ,RR(6 ,200) 
integer  error,  Hsubscript(33) 

common  /inv/Einv  /matrixE/E  /matrixA/A  /Hvector/H ,Hstar 
common  /matrixC/C  /loads/HH  /readings/RR  /misc/Rm 
data  Hsubscript/1 ,2,3,4,5,6,11,22,33,44,55,66,12,13,14,15,16, 
*  23,24,25,26,34,35,36,45,46,56,111,222,333,444,555,666/ 


c - 

C  Files  for  I/O:- 

C  SGBIN.DAT  =  file  containing  sets  of  H(i)  and  R(i)  for 

C  data  input 

C  SGB0UT.DAT  =  file  for  output  of  results 

C - 

open(unit  =  3,  f  ile=  ’  SGBIN .  DAT  ’  ,  status31  ’  old’ ) 
open(unit=4,  f ile= ’ SGBOUT .DAT ’  ,  status= ’new ’ ) 

C - 

C  Print  Headings  of  the  program  outputs 

C  - 

write(4 , 10) 

10  format (/  ,  ’  STRAIN  GAUGE  BALANCE  CALIBRATION  USING  LEAST’ 

*  ’  SQUARES  METHOD’/’  - ’ 

*  i - ’//’  The  Supplied  Component  Loads’ 

*  ’  and  Strain  Gauge  Bridge  Readings:-’/) 


M 


vrite(4 , 15)  (i,  i  =  l,6),  ( j ,  j  =  l,6) 

15  f ormat (6(5x , ’H’ ,il,3x) ,6(5x, ’R’ , i 1 ,3x)) 

0 - 

C  Initialize  the  [E]  and  [A]  matrices 

c - 

N  =  0 

do  n.roo  =  1,33 
do  n_col  =  1,33 

E(n_row ,n_col)  =  0.0 
end  do 
end  do 

do  n_col  =1,6 
do  n_row  =  1,33 

A(n_row ,n_col)  =0.0 
end  do 

sumR(i)  =  0.0 
end  do 

C  - 

C  Input  data  H(i)  =  SGB  readings  vector 
C  R(i)  =  Load  vector 

c  - 

20  read(3 , * ,end=50) 

read(3,»,end=50)  (H(i),i=l,6) 
read(3 , * ,end=50)  (R(i),i=l,6) 
write(4 ,30)  (H(i) ,i=l,6) ,  (R(i),i=l,6) 

30  format (12(x,f9.5)) 

N  =  N  +  1 
do  i  =  1,6 

sumR(i)  =  sumR(i)  +  R(i) 

c - 

C  Store  the  original  H(i)  and  R(i)  values  for  later  computations 

c - 

HH ( i , N )  =  H(i) 

RR(i ,N)  =  R(i) 
end  do 

C  - 

C  Evaluate  the  Least  Squares  elements  for  the  [E]  and  [A]  matrices 

C  - 

call  Hfactors(H) 
call  SumElements (R) 
go  to  20 
50  continue 

urite(4,60)  N 

60  format(’  Number  of  sets  of  loading  =  ’ 


,13,/) 


c  - 

C  Evaluate  the  inverse  of  [E] 

c - 

call  Matlnv(33,  E,  Einv,  error) 

if  (erro  .gt.  0)  go  to  80  !{  Ho  inverse  } 

C - 

C  Evaluate  the  coefficient  matrix  [C]  =  [Einv] [A] 

C - 

call  Mat Mult (Einv, 33 ,33, A, 33, 6, C, error) 
write(4,70)  (i,  i=l,6) 

70  f ormat ( 1H1 , / , ’  The  calibration  coefficient  matrix  [C] 

*4x ,  ’  H’ ’s’ ,6(7x, ’C’ ,il, ,4x)) 
do  i  =  1,33 

write(4,100)  Hsubscript ( i)  ,  (C(i,j),  j  =  l,6) 
end  do 
do  i  =  1,6 

Rm(i)  =  sumR( i)/f loat (N)  !{  mean  values  of  R(i)  } 

end  do 

call  Statistics (N) 
go  to  110 

80  write(* ,90) 

90  f ormat (/ ’  INVERSE  OF  [E]  DOES  NOT  EXIST  V.'/.V.’) 

100  f ormat (4x ,i3,6(2x,fl2.7)) 

110  close(4) 
end 


Hi 


Subroutine  SumElements (R) 

c  - 

C  Evaluate  the  sums  which  make  up  the  elements 
C  in  the  [E]  and  [A]  matrices 

c - 

double  precision  A(33,6),  E(33,33),  Hfactor(33) ,  R(6) 
common  /matrixE/E  /matrixA/A  /inv/Einv  /Hvector/Hf actor 
integer  row,  col 


do  row  =  1,33 
do  col  =  1,33 
E(row,col)  = 
end  do 
end  do 

do  col  =  1,6 
do  ion  =  1,33 
A(row,col)  = 


end  do 


E(row,col)  +  Hf actor (row) *Hf actor(col) 


A(row,col)  +  R(col) *Hfactor (row) 


end  do 
return 
end 


Subroutine  Hfactors(H) 

C  - 

C  Calculate  the  non-linear  load  factors 


double  precision  H(6) ,  Hfactor(33) 
common  /Hvector/Hf actor 

do  i  =  1,6 

Hfactor(i)  =  H(i)  !{ 

Hfactor(i+6)  =  H(i)*H(i)  !{ 

Hfactor(i+27)  =  H(i) *H( i) *H( i)  !{ 
end  do 
i  =  13 
do  j  =  1 , 5 

do  k  =  (  j  + 1 )  ,  6 

Hfactor(i)  =  H(j)+H(k)  !{ 

i  =  i  +  1 
end  do 
end  do 
return 
end 


linear  } 
squared  } 
cubed  } 


cross  product  } 


IX 


Subroutine  Statistics(N) 

C  - 

C  This  subroutine  evaluates  the  various  statistics  relating  to 
C  the  R(i)  values  as  estimated  by  the  calibration  equation 

C 

C  Input :  N  =  number  of  sets  of  loading 

C 

C  Associate  routines  SumStats(j) 

C  Statslnitialize 

C  SGBReading(Cl ,C2 ,Re ,H , error ) 

C  Tstats (ndf , t ) 

C  PartCorr(N ,i ,cpc) 

C  Chart ( scale ,N ) 

C  - 

double  precision  C(33 , 6 ) ,C1 (6 , 6) , C2 (6 , 27 ) ,R(6) ,H(6) ,Re (6) ,Rm(6) 
double  precision  HH(6 , 200) , RR(6 , 200) 

double  precision  sumRH(6 , 6) , sumR(6) ,sumH(6) , sumR2(6) , sumH2(6) 
double  precision  sumHH(7 , 7) ,Se (6) , SdR(6) ,Sdli (6) ,cmc (6) ,cl (6 ,6) 
double  precision  cpc (6 , 6) ,dRRm2 (6) ,dReRm2(6) , dRRe2(6) , res (6 ,200) 
common  /Sums/ sumRH , sumR , sumH , sumR2 , sumH2 , sumlffl , 

«  dRRe2 ,dReRm2 ,dRRm2 

common  /loads/HH  /readings/RR  /misc/Rm ,R ,Re ,H 
common  /residual/res .resmax  /matrixC/C 
parameter  (nv=33)  !{  no.  of  variables  } 

resmax  =  1.0e-20  !{  maximum  res  (residual)  for  chart  scaling  } 

call  Statslnitialize  !{  Initialize  the  variables  to  zero  } 

C  - 

C  Break  up  the  coefficient  matrix  into 
C  linear  (Cl)  and  non-linear  (C2)  matrices 
c  - 1 - 

do  i  =  1 ,  6 
do  j  =  1,6 

C 1 ( i  ,  j  )  =  C(j,i) 
end  do 
do  j  =  1,27 

C2 ( i  ,  j  )  =  C ( j  +6  ,  i ) 
end  do 
end  do 

C  -  - 

C  Restore  the  original  H(i)  and  R(i)  values  for  computations 

C  - 

do  j  1  ,  H 
do  i  -  1,6 

H(i)  -  HH(i ,j) 

R ( i )  ^  RR( i , j ) 
end  do 


c  - 

C  Estimate  the  values  of  R,  Re(i),  from  the  coefficients  matrices 
C  Cl  and  C2 ,  and  evaluate  the  various  sums  for  statistical  calculations 
C - 

call  SGBReading(Cl ,C2, Re, H, error) 
call  SumStats(j) 

C  - 

C  Use  HH(i,j)  as  storage  for  Re(i)  for  later  print-outs 

C  - 

do  i  =  1 ,  6 

Hll(i,j)  =  Re(  i ) 
end  do 
end  do 

20  ndf  =  float(N-nv)  !{  degree  of  freedom  } 

if  (N  . le .  nv)  ndf  =  1 
fN  -  fioat(K) 

call  Tstats (ndf , t)  !{  find  the  t-statistics  } 

c  - 

C  Evaluate  the  statistics 
C  Se(i)  =  standard  error  of  estimates 
C  SdR(i)  =  standard  deviation  of  R(i) 

C  SdU(i)  =  standard  deviation  of  H(i) 

C  cmc(i)  =  coefficient  of  multiple  correlation 
C  cpc(ij)  =  coefficient  of  partial  correlation 
C  cl(ij)  =  9S'/.  confidence  level  of  the  linear  coefficients 

C  - 

do  i  =  1,6 

Se(i)  =  sqrt (dRRe2 ( i ) /ndf ) 

SdR(i)  =  sqrt (dRRm2(i) /ndf ) 

SdH(i)  =  sqrt ( (sumH2(i)-sumH(i) *sumH( i)/f N)/ndf ) 
cmc(i)  =  sqrt (dReRm2{ i ) /(dRRe2{ i ) +dReRm2 ( i ) ) ) 
call  PartCorr(N  ,  i  ,cpc) 
end  do 
do  j  =  1,6 
do  i  =  1,6 

cl(i,j)  *  t/sqr t ( f loat (ndf ) ) *Se( i ) /SdH ( j ) 
end  do 
end  do 

C  - - - - 

C  Write  the  results  to  the  output  file 

C  -  -  •  - 

wr i te ( 4 , 30 ) 

write(4,110)  (i,  i *  1 , 6 ) 
write (4, 60)  (Se(i),  i=l,6) 

30  f  ormat  (  1 H 1  , /// ,  ’  Standard  error;-  of  •  •  :  rate'- :  ’,/) 

wr i te(4 ,40) 

wri  ted  ,110)  (  i  ,  1-1,6) 


write(4,60)  (cmc(i),  i=l,6) 

40  forraat(///,’  Coefficients  of  multiple  correlation:-’/) 
urite(4 ,50) 

write(4,120)  (i,  i=l,6) 

50  format(///,’  Coefficients  of  partial  correlation;-’/) 
do  j  =  1 ,  6 

write(4 ,65)  j,  (cpc(i,j),  i  =  1,6) 
end  do 

60  f ormat (6(2x  ,f 11 . 5) ) 

65  format(3x, ’H’ ,il ,6(2x,f 1 1 .5) ) 
write(4,70)  (i,  i=l,6) 

70  format(///  )5'/.  confidence  level  of  the  linear  coeff.:- 

*  2x,6  _>x ,  ’  C  ’  ,il ,  ’  -  ’  ,4x ) ) 
do  j  =  1 ,  6 

write(4,80)  (cl(i,j),  i  =  1,6) 
end  do 

80  format (6(2x,fl0.7)) 

C  - 

C  Print  out  the  measured  .estimated  and  the  residuals  of  R(i) 

c - 

write(4,90)  (i,  i=l,6) 

90  f ormat (1H1,’  Table  of  residuals:-  ’,//, 

*  6(8x , ’R’ ,il , llx)/6( ’  meas .  est .  res.  ’)) 
do  j  =  1 ,  Jf 

write(4,100)  (RR( i , j ) ,HH(i , j ) , res ( i , j )  ,  i  =  l,6) 

100  format (5(2(x,f5.3) ,x , f 6 . 4 , x , ’ I ’ ) , 2(x ,f 5 . 3) , x , f 6 . 4) 
end  do 

110  f ormat (3x ,6(6x , ’ R’ ,il ,5x) ) 

120  f ormat (8x ,6(6x , ’R’ , il ,5x) ) 

c  - 

C  Construct  the  chart  of  residuals 

c - 

scale  =  10.0/resmax 
call  Chart (scale ,N) 


’// 


return 

end 


Subroutine  Statslnitialize 

C - 

C  Initializes  the  variables  of  the  Statistics  routine  to  zero 

c - 

double  precision  sumRH(6 ,6) , sumR(6) , sumH(6) , sumR2(6 ) , sumH2 (6) 
double  precision  sumHH(7 , 7) ,dRRe2(6) , dRRm2 (6) tdReRm2 (6) 
common  / Sums/sumRH , sumR , sumH , sumR2 , sumH2 , sumHH , 

*  dRRe2 ,dReRm2 ,dRRm2 


do  i  =  1,  6 
do  j  =  1,6 

sumRH( j , i)  =  0.0 
sumHH( j +1 , i  +  1)  =  0.0 
end  do 

dRRe2(i)  =  0.0 
dReRm2(i)  =0.0 
dRRm2 ( i )  =  0.0 
sumR(i)  =  0.0 
sumH(i)  =  0.0 
sumR2(i)  =  0.0 
sumH2(i)  =  0.0 
end  do 
return 
end 


•)•> 


Subroutine  SGBReading(Cl , C2 ,R ,H, error) 

c - 

C  This  subroutine  evaluates  the  six  components  of  Strain  Gauge 
C  Balance  Readings  (R)  from  the  six  components  of  Loads  (H) . 

C 

C  Input:  C1,C2  -  linear  and  non-linear  coefficient  matrices 
C  Output:  R  (the  SGB  Reaings  vector) 

C  Error:  if  error  >  0,  an  error  has  occurred 
C 

C  Associate  routines:-  Hfactors(H) 

C  Mat Mult (Cl ,m,n,C2,n,k,C3, error) 

C - 

double  precision  Cl (6 , 6) ,C2(6 , 27 ) ,X(6) ,R(6) ,H(6) .Hstar (2T) 

integer  error 

common  /Hvector/X , Hstar 

call  Hf actors (H) 

C  - 

C  Perform  the  matrix  multipications : 

C  [Cl]  [H]  and  [C2]  [H*] 

C  - 

Call  MatMult(Cl ,6 ,6 ,H,6 , 1 ,R, error) 
if  (error  .gt.  0)  return 
Call  MatMult(C2,6,27 , Hstar, 27 , 1 ,X, error) 
if  (error  .gt.  0)  return 
do  i  =  1,6 

R(i)  =  R(i)  +  X(i) 
end  do 
return 
end 


Subroutine  SumStats(j) 


C - 

C  Evaluate  the  sums  for  the  statistical  estimations 


double  precision  sumRH(6 , 6) , sumR(6) , sumH(6) , sumR2 (6) , sumH2 (6) 
double  precision  sumHH(7 , 7) , res(6 , 200) , dRRe2 (6) , dRRm2(6) ,dReRm2(6) 
double  precision  H(6) ,R(6) ,Re(6) ,Rm(6) 
common  /Sums/sumRH , sumR , sumH , sumR2 , sumH2 , sumKH , 

*  dRRe2 ,dReRm2 ,dRRm2 

common  /misc/Rm,R,Re ,H  /residual/res , resmax 

do  i  =  1,6 

res(i,j)  =  R(i)-Re(i) 

if  (abs (res (i ,  j ) )  .gt.  resmax)  resmax  =  abs (res (i  ,  j )  ) 

dRRe2(i)  =  dRRe2(i)  +  res(i,j)**2 

dRRm2 ( i )  =  dRRm2(i)  +  (R(i)-Rm(i) )**2 

dReRm2(i)  =  dReRm2(i)  +  (Re(i) -Rm( i) ) **2 

sumR(i)  =  sumR(i)  +  R(i) 

sumH(i)  =  sumH(i)  +  H(i) 

sumR2(i)  =  sumR2(i)  +  R(i)*R(i) 

sumH2(i)  =  sumH2(i)  +  H(i)+H(i) 

do  m  =  1,6 

sumRH(i,m)  =  sumRH(i,m)  +  R(i)*H(m) 
sumHH(i+l ,m+l)  =  sumHH(i+l ,m+l)  +  H(i)*H(m) 
end  do 
end  do 
return 
end 


c 


Subroutine  Tstats (ndf , t ) 


C  Estimate  the  t-statistics  from  the  degree  of  freedom 

c  - 

dimension  tdf(30) 

data  tdf/12. 706, 4. 303, 3. 182, 2. 776, 2. 571, 2. 447, 2. 365, 2. 306, 

*  2.262,2.228,2.201,2.179,2.160,2.145,2.131,2.120, 

+  2.110,2.101,2.093,2.086,2.080,2.074,2.069,2.064, 

*  2.060,2.056,2.052,2.048,2.045,2.042/ 


if  (ndf  .le.  30)  then 
t  =  tdf(ndf) 

else  if  (ndf  . le .  40)  then 

t  =  tdf (30) -0.021 *f loat ((ndf -30)/ 10) 
else  if  (ndf  .le.  60)  then 

t  =  2.021-0. 021*f loat ( (ndf -40) / 10) 
else  if  (ndf  . le .  120)  then 

t  =  2.000-0. 020*f loat ( (ndf -60)/ 10) 
else 

t  =  1.96 
end  if 
return 
end 


Subroutine  PartCorr ( N , iR , cpc ) 


C - 

C  Subroutine  to  evaluate  the  coefficients  of  partial  correlation 
C 

C  Inputs:-  N  =  total  number  of  sets  of  data 

C  iR  =  the  subscript  of  R  chose  coeff.  of  part.  corr. 

C  are  to  be  evaluated 

C  Outputs:-  cpc  =  arrays  containing  the  coeff.  of  part.  corr. 

C 

C  Cautions:-  Errors  may  occur  if,  during  the  evaluation  process, 

C  the  values  of  certain  correlation  coefficients  are 

C  close  to  1 . 

C 

C  Associated  routine  :-  Rahk(m,n) 

C - 

double  precision  sumRH(6 , 6) , sumR(6) , sumH(6) , sumR2(6) , sumH2 (6) 
double  precision  sumXX (7 , 7)  ,  sumX (7)  ,  sumX2 (7) 
double  precision  corr (6 ,6 ,7) , cpc(6 ,6) , a ,b , c 
integer  mediator(6) ,  order 

common  /Sums/sumRH , sumR , sumH , sumR2 , sumH2 , sumXX 


c - 

C  Map  the  sums  of  R(iR)  and  H(i)  into  the  array  sumX(j) 

C  the  sums  of  squares  into  the  array  sumX2(j),  and 

C  the  sums  of  cross  products  into  sumXX( j ,k) 

c  - 

sumX(l)  =  sumR(iR) 
sumX2(l)  =  sumR2(iR) 
do  j  =  2,7 

sumX(j)  =  sumH(j-l) 
sumX2(j)  =  sumH2(j-l) 
sumXX(l,j)  =  sumRH(iR , j-1 ) 
end  do 

c  - 

C  Evaluate  the  zero  order  (order  =  1)  correlation  coefficients  first 

c - 

do  ml  =  1,6 

do  m2  =  (ml  + 1 )  , 7 

corr ( 1 ,ml ,m2)  =  ( N*sumXX (ml ,m2) -sumX (ml ) +sumX (m2) ) 

*  /sqrt((H*sumX2(ml ) -sumX(ml ) *sumX (ml ) ) * 

*  (N*sumX2(m2)-sumX(m2) *sumX(m2) ) ) 
end  do 

end  do 


2(5 


C - 

C  Start  evaluating  higher  order  correlation  coefficients 

C - 

do  i  =  2,7  !{i=  component  subscript  } 

call  rank(i .mediator) 

do  order  =  2,6  !{  order  of  correlation  coefficient  } 

med  =  mediator (order )  !{  subscript  of  the  mediator  } 

do  nl  =  l.(med-l)  !{  subscript  of  the  1st  variable  } 

do  n2  =  (nl+l),7  !{  subscript  of  the  2nd  variable  } 

if  ((n2  .It.  med)  .or.  (n2  .eq.  i))  then 
a  =  corr(order-l ,nl ,n2) 
b  =  corr (order-1 ,nl , med) 
c  =  corr(order-l,n2,med) 

if  (n2  .gt.  med)  c  =  corr(order-l ,med,n2) 

c - 

C  To  prevent  negative  sqrt  argument  and  also  divide  by  zero  errors 
C  the  correlation  coefficient  is  set  slightly  less  than  1  if 
C  either  b  or  c  >=  1. 

C - 

if  (b  .ge.  1.0)  b  =  0.999999 

if  (c  .ge.  1.0)  c  =  0.999999 

corr (order ,nl ,n2)  =  (a-b*c)/sqrt ( ( 1 .0-b*b) * ( 1 . 0-c*c) ) 
end  if 

end  do  ! {  end  of  n2  loop  } 

end  do  !{  end  of  nl  loop  } 

end  do  ! {  end  of  order  loop  } 

end  do  ! {  end  of  i  loop  } 

do  k  =  1,6 

cpc(iR.k)  =  corr(6 , 1  ,k+l) 
end  do 
return 
end 
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Subroutine  Rank(n,med) 


C - 

C  Rank  the  subscripts  of  the  mediator  from  7  down  to  2  excluding  n 

C - 

integer  med(6) 

med(2)  =  7 

if  (n  .eq.  7)  med(2)  =  med(2)-l 
do  i  =  3,6 

med(i)  =  med(i-l)  -  1 

if  (med(i)  .eq.  n)  med(i)  =  med(i)  -  1 
end  do 
return 
end 


28 


Subroutine  Chart(scale ,N) 


c - 

C  Construct  and  print  the  chart  of  the  residuals  of  the  difference 
C  between  the  estimated  and  measured  values  of  R(i) 

C 

C  Inputs:-  scale  =  10/resmax 

C  N  =  no.  of  sets  of  data 

C  Outputs:-  line  printer  plot  of  the  chart  of  residuals  written 
C  to  logical  unit  number  4 

C 

C  Associate  routine  :-  Resetline(line) 

c - 

double  precision  res (6, 200) 
character  line(131) ,  symbol(6) 
data  symbol/ ’  *  ’  ,  ’ o’, o’ 
common  /residual/res , resmax 

write(4,5)  resmax,  (i,  i=l,6) 

5  format(lHl,’  Deviations  of  the  measured  R’’s  from’ 

*  ’  the  estimated  R’ ’s’/, ’  Maximum  deviation  =  ’,el2.5,//, 

*  6(l0x,’R’,il,10x)) 
do  i  =  1 ,  M 

call  resetline(line)  !{  clear  the  line  of  characters  } 
do  j  =  1,6 

index  =  res( j ,i) *scale 
index  =  index  +  (j-l)*22  +  11 
line(index)  =  symbol(j) 
end  do 

write(4,15)  (line(k) ,  k=l,131) 

15  format(131al) 
end  do 
return 
end 
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Subroutine  Resetline(line) 


C - 

C  Clears  a  line  of  character  to  contain  spaces  and 
C  the  zero  marker  ’ ! 1  only 
C  -  Called  by  the  subroutine  Chart 

C - 

character  line(13l) 

do  i=l,  131 
lined)  =  1  ’ 
end  do 

do  i  =  11,131,22 
lined)  =  ’  !  ’ 
end  do 
return 
end 


:sn 


Subroutine  Mat Mult (ml , rowl,coll,m2,row2,col2 .product .error) 

C - 1 

C  Matrix  multipication  routine  ! 

C  ! 

C  Input:  ml,  m2  (the  first  and  second  matrices)  ! 

C  rowl,  coll,  row2,  col2  (dimensions  of  the  two  matrices)  ! 

C  Output:  product,  error  ! 

C  ! 

C  Error:  0  -  successful  operation  ! 

C  1  -  the  column  of  the  first  matrix  does  not  match  that  of  ! 

C  the  second  ! 

c - . 

integer  rowl,  row2,  coll,  col2,  error 

double  precision  ml (rowl , coll ) ,m2(row2 ,col2) .product (rowl , col2) 
if  (coll  -  row2)  10,  20,  10 

c - 

C  No.  of  columns  in  ml  does  not  match  that  in  m2  -  no  operation 

C - 

10  error  =  1 
return 

c - 

C  Perform  the  Matrix  Multipications 

C  - 

20  do  i  =  1 ,  rowl 

do  j  =  1,  col2 

product(i.j)  =  0.0 
do  k  =  1,  coll 

product(i.j)  =  product(i.j)  +  ml (i ,k) *m2(k , j ) 
end  do 
end  do 
end  do 
error  =  0 
return 
end 


:si 


Subroutine  Matlnv(dimen,  matrix,  inverse,  error) 


c  — 

— 

c 

Purpose  : 

:  To 

calculate  the  inverse  of 

a  matrix 

c 

with  Gauss-Jordan  elimination. 

c 

c 

Input  : 

:  dimen,  matrix 

c 

Output 

:  inverse,  error 

c 

c 

Remarks  : 

:  The  original  matrix  will  be 

changed  ' 

c 

an 

identity  matrix. 

c 

c 

Error  code 

:  0 

-  no  error 

c 

1 

-  dimension  of  matrix  <  1 

c 

2 

-  no  inverse  exists 

c 

c 

Written  by 

:  S 

.  Lam 

c 

Date 

:  2 

March  1988 

c - • 

integer  dimen,  error,  row,  ref_row,  term 

double  precision  matrix(dimen , dimen) , inverse (dimen, dimen) 
double  precision  almostzero/1 . 0e-20/ ,  multiplier 
common  /zero/almostzero 


call  Initial (dimen,  matrix,  inverse,  error) 
if  (dimen  .le.  1)  return 
ref_row  =  0 

do  while  ((error  .eq.  0)  .and.  (ref_row  .It.  dimen)) 
ref_row  =  ref_row  +  1 


Check  to  see  if 

the  diagonal  element  is  zero 

+ 

if  (abs(matrix(ref _row,ref _row) )  .It.  almostzero) 
call  Pivot(dimen,  ref_row,  matrix,  inverse,  error) 
if  (error)  30,  20,  30 

If 

error  =  0  , 

carry  out  the  inverse  operation 

20 

divisor  = 

matrix(ref_row,refrow) 

do  term  =  1 ,  dimen 

matrix(ref _row  ,  term)  =  matrix ( ref _row , term) /divisor 
inverse(ref _row , term)  =  inverse (ref _row , term) /divisor 
end  do 


:v> 


do  row  =  1 ,  dimen 


C 

C 

C 


Make  the  ref_row  element  of  this  row  zero 


if  ((row  .ne.  ref_row)  .and. 

(abs (matrix(row  ,  ref _row) )  ,gt.  a1 mostzero) )  then 
multiplier  =  -matrix(row ,ref_row)/matrix(ref _row ,ref _row) 
do  term  =  1 ,  dimen 

matrix(row .term)  =  matrix (row , term)  +  multiplier* 
matrix ( ref _row , term) 

inverse (row , term)  =  inverse (row , term)  +  multiplier* 
inverse (ref _row .term) 

end  do 
end  if 
end  do 


C  If  error  <>  0,  terminate  operation 

C - 

30  end  do  !{  end  of  do  while  loop  } 
return 


end 


:n 


Subroutine  Ini t i al (dimen ,  matrix,  inverse,  error) 

C  - ! 

C  This  subroutine  tests  for  errors  in  the  value  of  dimen  ! 

C  ! 

C  Input:  dimen,  matrix  ! 

C  Output:  inverse,  error  ! 

c  - j 

integer  dimen,  error,  row,  col 

double  precision  matrix(dimen, dimen) ,  inverse(dimen, dimen) 
double  precision  almostzero 
common  /zero/almostzero 

error  =  0 

if  (dimen  -  1)  10,  20,  30 

C  - 

C  Case  dimen  <  1 ,  error 

0 - 

10  error  =  1 
return 

c  - 

C  Case  dimen  =  1,  single  element  matrix 

c  - 

20  continue 

if  (abs (matrix( 1 , 1) )  .It.  almostzero)  then 
error  =  2  {  singular  matrix  } 

else 

inverse(l.l)  =  1 . 0/matrix ( 1 , 1 ) 
end  if 
return 

c - 

C  Case  dimen  >  1,  make  the  inverse-to-be  an  identity  matrix 

c  - 

30  do  col  =  1,  dimen 

do  row  =  1 ,  dimen 

if  (col  .eq.  row)  then 


inverse (row , 

col ) 

=  1.0 

else 

inverse ( row , 

col ) 

=  0.0 

end  if 
end  do 


end  do 
return 
end 


Subroutine  Pivot(diraen,  ref_row,  matrix,  inverse,  error) 

c - 

C  This  subroutine  searches  the  ref_row  column  of  the  given  matrix 

C  for  the  first  non-zero  element  below  the  diagonal.  If  it  finds 

C  one,  then  the  subroutine  switches  rows  so  that  the  non-zero  element 
C  is  on  the  diagonal.  This  same  operation  is  applied  to  the  inverse 
C  matrix.  If  no  non-zero  element  exists  in  a  column,  the  matrix  is 
C  singula*-  and  no  inverse  exists. 

C 

C  Input:  dimen,  ref_row,  matrix,  inverse 
C  Output  :  matrix,  inverse,  error 

c  - - 

integer  dimen,  error,  neo.row 

double  precision  matrix(dimen, dimen) ,  inverse(dimen, dimen) 
double  precision  almostzero 
common  /zero/almostzero 

error  =2  !{  No  inverse  exists  } 

new_row  =  ref_row 

c  - 

C  Try  to  find  a  row  with  a  non-zero  diagonal  element 

C  - 

do  while  ((error  .gt.  0)  .and.  (new_row  .It.  dimen)) 
new_row  =  new_roa  +  1 

if  (abs(matrix(new_row,ref_row))  .gt.  almostzero)  then 

C  - 

C  Switch  rows  between  the  new_row  and  the  ref_row 

c  - 

call  SwitchRow (matrix ,  dimen,  new_row,  ref_row) 
call  SwitchRow(inverse ,  dimen,  new_row,  ref_row) 
error  =  0 
end  if 

end  do  !{  While  } 

return 

end 


:sr> 


Subroutine  Sw itchRow (matrix ,  dimen,  roul  ,  row2) 


c  - 

C  This  subroutine  switches  the  rows  specified  by  rowl  and  row2 
C  of  the  given  matrix. 

C 

C  Input:  matrix,  dimen,  rowl,  row2 

C  Output:  matrix 

c  - 

integer  dimen,  rowl,  row2,  term 

double  precision  matrix(dimen , dimen) ,  dummy 


do  term  =  1 ,  dimen 
dummy  =  matrix(rowl 
matrix(rowl .term)  = 
matrix(row2 , term)  = 
end  do 
return 
end 


.term) 

matrix(row2 

dummy 


term) 


6  Example 


As  a  demonstration  of  the  functions  of  the  program  SGBCalib,  a  “numerical”  calibration 
was  carried  out  using  an  “artificial”  balance  whose  calibration  coefficient  matrix  is  given  as 
follows: 


This  represents  a  perfectly  linear  balance  with  no  interactions  between  the  components.  It  is 
an  ideal  balance  which,  of  course,  is  practically  impossible  to  attain.  Nevertheless,  it  serves 
the  purpose  of  a  simple  illustration  of  the  capabili  ies  of  the  calibration  program. 

:(7 


A  number  of  sets  of  randomly  generated  component  loads  //,,  (ranged  from  0  to  2 
units)  were  used,  in  conjunction  with  the  coefficient  matrix  (C]  of  equation  (181,  to  produce  a 
number  of  corresponding  sets  of  strain  gauge  balance  outputs  /<!,  according  to  the  relationship 
described  by  equation  (2).  The  resulting  /f,  values  were  multiplied  by  a  random  factor  in 
the  range  of  ±0.002  to  simulate  a  random  noise  level  of  ±0.2%  on  the  strain  gauge  balance 
outputs.  The  sets  of  Hi  and  noise  inflicted  /2,  values  were  written  to  the  input  file  SGBIN.DAT 
according  to  the  format  described  in  Section  4.2.1.  The  program  was  then  executed  to 
evaluate  the  elements  of  [C]  from  these  sets  of  If ,  and  /£,-  values.  The  outputs  of  the  program 
are  shown  in  page  39  through  to  page  43. 

Examination  of  the  results  shows  that  the  maximum  deviation  of  the  estimated  calibra¬ 
tion  coefficients  is  approximately  1%  from  the  “true”  values.  This  is  hardly  acceptable  in 
the  high  precision  work  that  a  strain  gauge  balance  is  expected  to  operate.  The  inaccurate 
estimation  of  the  coefficients  in  this  particular  case  is  because  not  all  of  the  randomly  applied 
component  loads  necessarily  give  useful  independent  information  for  the  evaluation  of  the 
coefficients.  The  ±0.2%  random  noise  level  further  distorts  this  information  from  accurately 
determining  the  coefficients.  The  accuracy,  however,  could  be  improved  if  a  larger  number  of 
sets  of  loadings  is  used.  Pages  44  and  45  are  an  abstract  of  the  outputs  from  a  similar  test  run 
except  that  100  sets  of  random  loadings  were  used  for  the  “calibration”.  The  results  showed 
considerable  improvement  in  the  accuracies  of  estimation  of  the  calibration  coefficients. 

On  the  other  hand,  the  correlation  coefficients  show  perfect  multiple  correlations  between 
the  strain  gauge  balance  outputs  and  the  six  component  loads.  There  is  also  perfect  linear 
relationship  between  each  of  the  R{  and  its  corresponding  H{.  These  are  to  be  expected 
because  the  R{  values  were  derived  from  exactly  the  same  relationship  as  the  regression 
equation. 
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